// ANSWER: 269533451410884183

#include <iostream>
#include <algorithm>
#include "prime.h"

// Computes the sum of f(k^3) for 1 <= k <= K.
static long long solve(int K, std::vector<long long> &terms)
{
	// Let u[k] be the greatest prime factor of (k+1).
	std::vector<int> u(K+1);
	for (int k = 1; k <= K; k++)
	{
		if (u[k] == 0) // is prime
		{
			int p = k + 1;
			for (int j = k; j <= K; j += p)
				u[j] = p;
		}
	}

	// Let v[k] be the updated version of k(k-1)+1.
	// Let largest_prime[k] be the largest prime divisor of v[k].
	std::vector<long long> v(K+1);
	std::vector<long long> largest_prime(K+1);
	for (int k = 1; k <= K; k++)
	{
		v[k] = (long long)k*(k-1)+1;
	}
	largest_prime[1] = 1;
	for (int k = 1; k <= K; k++)
	{
		long long p = v[k];
		if (p > 1)
		{
			for (long long j = k; j <= K; j += p)
			{
				if (largest_prime[j] < p)
					largest_prime[j] = p;
				do
				{
					v[j] /= p;
				}
				while (v[j] % p == 0);
			}
		}
	}

	// Now find the sum of f(k^3).
	long long sum = 0;
	for (int k = 1; k <= K; k++)
	{
		long long p = std::max((long long)u[k], largest_prime[k]);
		sum += (p - 1);
		terms[k] = (p - 1);
	}
	return sum;
}

// Computes the sum of f(k^3) for 1 <= k <= K in a brute force way.
static long long verify(int K, std::vector<long long> &terms)
{
	std::vector<int> primes;
	find_primes_below(K+2, primes);

	long long sum = 0;
	for (long long k = 1; k <= K; k++)
	{
		long long m = 1;
#if 1
		prime_factorize(k*k*k+1, [&m](long long p) { m = p; }, primes.cbegin(), primes.cend());
#else
		prime_factorize(k*k*k+1, [&m](long long p) { m = p; });
#endif
		sum += (m - 1);
		terms[k] = (m - 1);
	}
	return sum;
}

static void show_factors(long long n)
{
	std::cout << n << " =";
	prime_factorize(n, [](long long p) { std::cout << " " << p; });
	std::cout << std::endl;
}

void solve_problem_343()
{
#if 0
	const long long K = 10000;
#else
	const long long K = 2000000;
#endif

	std::vector<long long> terms(K+1);
	std::cout << solve(K, terms) << std::endl;

#if 0
	std::vector<long long> terms2(K+1);
	std::cout << verify(K, terms2) << " - verify" << std::endl;
	for (int k = 1; k <= K; k++)
	{
		if (terms[k] != terms2[k])
		{
			std::cout << "Mismatch: " << k << " " << terms[k] << " " << terms2[k] << std::endl;
			show_factors((long long)k*k*k+1);
		}
	}
#endif
}
